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I. INTRODUCTION 


Two equation turbulence models used in the computation of the turbulent flows involve the 
specification of five or six constants. These constants have been evaluated using the experimental data 
on boundary layer flows. These are not changed during any calculation. However, these are not 
universal constants. They need to be changed for predictions of different classes of flows, for example 
predictions of plane two-dimensional jet and axisymme trie jet [1,2] The constants have to be changed 
to correctly predict the spread rate of an axisymmetriejet. The constants have to be changed in order 
to accommodate the effects such as streamline curvature, low Reynolds number, near wall region, etc. 
[31 


Another important element in the predictions is the specification of inlet boundary conditions. 
Sometimes they are considered of minor importance while evaluatmg the performance of turbulence 
models against experimental data. Assumed boundary conditions can lead to wrong conclusions about 
the performance of the model [3,4]. In particular, the inlet profiles of k and e can have a significant 
effect on the flow downstream. The profiles of k may be evaluated from the experimental data. But 
no satisfactory method of evaluating the inlet profile of e is available The inlet conditions become 
more controlling when the inlet flow has a swirl component. For example, only with a “suitable” initial 
condition, can the coswirl flow in coaxial jets be predicted [5] In addition, minor changes in the e 
level of the inlet section have been found to have the same degree of modifications on the flow field as 
the swirl related corrections to the turbulence model [6] 

In the case of complex flow configurations, such as the SSME internal flow configuration, it is 
not expected that detailed inlet boundary conditions will be available from experimental measurements. 
In such cases it is important to include a sensitivity test to the initial conditions in turbulence model 
predictions. The issue of initial conditions is more important in internal flows, where elliptic procedures 
are used for the numerical solutions, than the free shear flows or boundary layer type flows, where 
parabolic, marching techniques are employed. In boundary layer type flows, the evaluation of the 
physical turbulence model is always made in the self preserving region, where the influence of the inlet 
condition is insignificant [7] In internal flows, especially flows with recirculation, the influence of the 
inlet boundary conditions may extend up to the reattachment point. 

In this report an effort is made to understand the difference in k-e model predictions of confined 
plane two-dimensional expansion flow, and axisymmetric sudden expansion flow. Model constants are 
not changed to correctly predict the experimental data. Only the predictions for the same initial flow 
conditions and expansion ratio are examined. The predicted results of the flow over a backward-lacing 
step, the symmetric plane two-dimensional expansion flow, and the axisymmetric expansion flow are 
compared. The report also presents the calculations of the flow in a confined coaxialjet. The predic- 
tions of the coaxialjet with the velocity ratio of annulus to central jet, equal to unity are compared 
with those of an axisymmetric sudden expansion flow. Results are also presented for different ratios 
of annulus to centraljet velocities. 



The effects of inlet k-e profiles and Reynolds number on the predicated reattached lengths in the 
case of axisymmetric expansion are investigated, The effects of inlet turbulence level are significant for 
the case of confined coaxial jet in which the development of the flow is highly dominated by the coupled 
diffusion processes and the interactions of the incoming annular and central jet flows. The situation is 
more complicated when the inlet swirl is present in the annular region where the diffusion process is 
coupled with the extra strain induced by the streamline curvature. An effort is also made to investigate 
the sensitivity of the inlet turbulence level on the swirling flow field. 

The predictions reported here have been earned out using PHOENICS [8] and TEACH [9] codes. 
These codes use in principle the same form of the equations governing the flow [10] However, they 
differ in the details of the solution procedures. A comparison of the prediction of the two codes for the 
axisymmetric sudden expansion flow is presented. 


II. PLANE TWO-DIMENSIONAL AND AXISYMMETRIC RECIRCULATING FLOWS 


The k-e turbulence model is used in a majority of all turbulent flow predictions, The standard 
model uses five constants determined for boundary layer flows. These constants are not as universal as 
one expects them to be. For example, the standard model predicts the spread rate of a plane jet 
correctly while it overpredicts the spread rate of a circular jet, Several modifications of the model have 
been proposed for improving the prediction of circular jets [1], However, no modification of the model 
produces correct predictions for all axisymmetric and plane two-dimensional cases [2]. 

A close examination of the confined recirculating flows seems to indicate that such anomaly 
between the plane two-dimensional and axisymmetric predictions is present in the computations of 
internal flows as well. Specifically, the length of the recirculation region, X r is underpredicted by about 

20 percent for the flow over a backward-facing step, whereas X f for an axisymmetric sudden expansion 

is predicted within the expenmental accuracy. The reattachment lengths reported by several investigators 
for backward-facing step and axisymmetric expansion, respectively, illustrate the point [31 The poor 
performance of the k-e model for the backward-facing step flow is not clear, 

A specific comparison of the prediction of the recirculation region for the two cases with 
measured data is shown in Figure 1. In Figure la, the length of the recirculation region is predicted 
accurately for the pipe expansion [ 12,131 The figure shows the locus of the flow reversal (U = 0). 

The flow reversal line shown for the backward-facing step (Fig. lb) shows the rather severe underpre- 
diction of the separation length [11] In both cases the width of the recirculation region is not very well 
predicted. 


A. Boundary Condition Effects 

In this report, an effort is made to explain the differences in predictions of plane and axisym- 
metric-confined recirculating flows. The flow over a backward-facing step, the flow in a plane symmetric 
expansion, and the flow in an axisymmetric pipe expansion are computed for the same inlet conditions. 

The Reynolds number based on the upstream channel width (diameter) is 9.3 x 10^ A grid of 42 x 42 
is used to produce grid independent solutions. The expansion ratio is 1.5 for which experimental data 
for backward-facing step flow is available. 
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(a) Axisymmetric expansion. 
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Figure 1 , Locus of flow reversal. 


1 The Keattachment Length 

The single most important parameter to compare in these flows is the reattachment length X r 

The reattachment length for the backward-facing step is obtained as 5.8H as compared to the “accepted” 
experimental value of 7H [11], where H is the step height. The X r for the pipe expansion is 7 7H for 

this expansion ratio. The symemtric two-dimensional expansion has an intermediate value of 6.3H, 

The reduction of X f from 6.3H for symmetric expansion to 5.8H for the asymmetric expansion is due 

to the stabilizing effect of the presence of the top wall instead of a symmetry line. The effect is not 
very significant in this case but has been found to be significant in the case of bluff body stabilized 
flows [14]. Now, consider only the plane two-dimensional symmetric expansion and axisymmetric 
expansion. In the limiting case of walls at infinity, free shear flow, the spread rate of the plane- 
separated layer is higher than the axisymmetric one, as shown by the experiments. For the finite step 
height the initial shear layer is unaware of the wall and thus the spread rates follow more or less the 
same trend as in the unconfined flows. Thus, for the plane symmetric expansion, X r is shorter than that 
for the axisymmetric expansion. 
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a Mean Velocity Profiles 

Prediction of the flow over a backward-lacing step was studied extensively in the 1980-81 
Stanford conference on complex turbulent flows. The mean velocity profiles behind the step at three 
different downstream distances obtained in the present study are shown in Figure 2. The experimental 
data for the backward-facing step and the predictions for plane symmetric expansion and the axisym- 
metric expansion are also shown in the figure. It is seen that the mean velocity profile in the separated 
region is predicted inaccurately. This seems to stem from the wall function used. Because of the nature 
of implementation of the wall function used, the first point next to the wall has the hghest return flow 
velocity, in all the computations. This is not true in the actual flow. The way of handling the wall 
function in the PHOENICS code appears to result in a peaked mean velocity from the near-wall node 
point to the subsequent points [15] This is discussed in Section IV 




X/H- 8 



Figure 2. Mean velocity profiles (X/H = 1,33 ? 5.33, 8.0). 


The mean velocity profile predictions for the plane symmetric expansion and the flow over a 
backward-lacing step show more or less a similar trend. In the axisymmetric expansion case, the flow 
(development) expansion is slower as shown by the mean velocity profiles (Fig. 2a-2c), This is also seen 
in the variation of the centerline velocity in the streamwise direction (Fig. 3), In the region just down- 
stream of the expansion, the decay rate of the centerline velocity is substantially higher for plane 
symmetric expansion, for a higher expansion ratio of 2 (Fig. 3b). The experimental data of Reference 12 
for the pipe expansion are also shown. The prediction is in fairly good agreement with the measure- 
ments. The faster decay of the centerline velocity in the plane flow causes a stronger streamline curva- 
ture than with the axisymmetric case. The computation of the production term in the e-equation 


9Uj 

p =~ u i u i . 

n 3x j 

tends to be more inaccurate in the plane expansion, For this reason, even the correction of the produc- 
tion and dissipation terms in the e-equation for the k-e model does not improve its performance for the 
plane flow The algebraic stress model employing an algebraic equation for ujuj, can only improve the 
prediction of the two-dimensional expansion [15,30], 



Figure 3. Variation of centerline velocity. 
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b< Kinetic Energy 

The kinetic energy profiles just downstream of the expansion are shown for the three cases in 
Figure 4. For all the three cases, the profiles show the characteristics of the separated region. That is, 
just downstream of the expansion, the kinetic energy maximum occurs in the separated layer at a y/H 
of 1.0. Further downstream of the expansion, the location of the kinetic energy maximum moves closer 
to the wall, as in the attached flows. However, the kinetic energy level is higher for the axisymmetric 
expansion than for the plane expansion flow, The kinetic energy contours in the entire computational 
field are shown in Figure 5 for the three cases, for an expansion ratio of 1.5. Figure 6 shows the con- 
tours for an expansion ratio of 2. For this expansion ratio the contours for the backward-facing step 
and the symmetric expansion are similar, except that the symmetric case shows a slightly shorter core 
region. The contours for the axisymmetric expansion clearly show a shorter core region as compared 
to the plane flow This is in accordance with the experimental observations [17] , Also, as observed 
in the experiments, the core region extends beyond the reattachment point in the plane flow, while in 
the axisymmetric flow the opposite holds (Fig. 7). 

X/H = 1.33 X/H = 5.33 



k/U 0 2 x 10 2 

Figure 4. Kinetic energy profiles (X/H = 1.33, 5.33, 8.0). 
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(a) Backward-Facing step. 



(b) Plane symmetric expansion. 




(c) Axi symmetric expansion. 

Figure 5 Contours of kinetic energy (expansion ratio = 1,5). 
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(b) Plane symmetric expansion. 



(c) Axi symmetric expansion. 

Figure 6. Contours of kinetic energy (expansion ratio = 2.0). 


c. Dissipation Rate of Kinetic Energy 

The kinetic energy dissipation rate profiles for an expansion ratio of 2 are shown in Figure 8. As 
for the kinetic energy, the dissipation rate contours for the backward-facing step and the symmetric plane 
expansion are similar - . However, the contours for the axisymmetric case show the short-core region. 

d. Wall Static Pressure Variation 

The wall static pressure variations for the three cases for an expansion ratio of 1 5 are shown in 
Figure 9 The pressure recovery downstream of the reattachment is faster in the case of plane two- 
dimensional expansions compared to the axisymmetric case. This is consistent with the shorter reattach- 
ment lengths predicted for plane flow 
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Figure 9. Wall static pressure distributions. 


B. Sensitivity of the Axisymmetric Expansion Flow to Inlet Boundary Conditions 

The inlet mean velocity profile is normally available from the measurements. Otherwise, either a 
fully developed flow or a plug flow is specified in the predictions. Measurements show that a thicker 
inlet boundary layer produces a larger recirculation region [30] . However, the predictions show only a 
small increase in X f (about 2 percent) for the fully developed flow as compared to plug flow 


The specification of inlet profiles for k and e is more difficult. The k profile may be estimated 
from the measurements. In general, k is specified as a percentage of the inlet mean square velocity 


k : 


in 


0.003 x U in 2 


( 1 ) 


Depending on the percentage (constant), the predicted results vary. For example the reattachment 
length X can vary as much as 1.0H as shown in Figure 10, by varying the constant, 
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io -4 io -3 io -2 io -1 io° 
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i / H or k/Uin 


Figure 10. Variation of X f with inlet k and e levels. 

The inlet profile for e has to be assumed, since no measurements are possible. It is specified as 

e = k^^/(Const £) (2) 

where £ is the length scale and = 0.09 The variation of the constant from 1.0 to 0.003 does not 

change the predicted reattachment length appreciably Only when the constant is specified as 0.0003, 
i.e., Q = 0.0003H, the reattachment length increases by about 0.8H. One point to note here is that the 
results show the correct expected trend. That is, with a larger inlet length scale (higher viscosity), the 
reattachment length is shorter For the case of the flow over a backward-facing step, such a trend was 
not observed [ 161 
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Moon and Rudinger [13] found experimentally that for Reynolds numbers greater than lCr , 
the length of the recirculation was independent of the Reynolds number, Red (based on the upstream 

mean velocity and diameter). The predicted results show that the reattachment length increased slightly 

with increase in Reynolds number at low Reynolds numbers (Fig. 11), and for Re^ > 10^ , no signific- 
and change in X r is observed as in the experimental results. 



Figure 11. Variation of X r with Reynolds number 


The effect of expansion ratio D/d on the length of the separated region is shown in Figure 12 
These results were obtained for a Reynolds number of 2.1 x 10^ The recirculation length X r is seen 

to increase with increase in expansion ratio for all the three cases. However, the variation is not linear 
as in the laminar flows. The difference in the values of Xj. for axisymmetric expansion and plane 



EXPANSION RATIO, D/d 

Figure 12. Variation of X r with expansion ratio. 
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symmetric expansion is highest (2.3H) for an expansion ratio of 1.5 At the expansion ratio of 2, the 
plane symmetric expansion will be unstable and unequal regions of separation on either side of the 
centerline occurs. Here, the solution has been obtained by considering only one -half of the duct. As the 
expansion ratio becomes larger (=2), the stabilizing effect of the top wall becomes significant. 


III. CONFINED FLOWS WITH/WITHOUT SWIRL 


The study of axisymmetric confined flows actually involves three cases. (1) pipe expansion, 

(2) coaxial jet in confmed sudden expansion without swirl, and (3) coaxial jet in confined sudden expan- 
sion with swirl in the outer annular jet. The flow in a pipe expansion can be viewed as a coaxial jet 
having equal velocities between the annularjet and centraljet, with the dividing lip of the central pipe 
being infinitely thin. Table 1 summarizes some of the important parameters of the cases studied. In 
these cases, the ratio of the velocity of the annularjet, U a> to the velocity of the centraljet velocity, 

Up increases from 1 to 12. The difference between the central and annularjet velocities in a confined 

chamber will create different turbulence diffusion processes. The standard k-e model is used for the 
predictions of these flows with the same inlet flow turbulence level as calculated from equations (1) 
and (2). The sensitivity study of the initial boundary conditions on the flow fields for a specific case 
(U a /Up = 3) will be presented in Section 1II.B, 

TABLE 1 CONFINED COAXIAL JET WITH DIFFERENT U a /U p 


Case 

Reynolds Number 


*V D P 

D o/ D a 

Exp. X r /H 

Predicted X r /H 

Moon and Rudinger [13] 

2.611 x 10 5 

1 


1.43 

8.7 

7.2 


Re = 39900 
Re P = 55200 

1 



9.25 

8.2 

Habib and Whitelaw [19] 

Re p = 20600 
Re, = 85000 

d 

3 

2.76 

2.81 

10.25 

8.4 

Johnson and Bennett [20] 

Re = 15900 
Re P =47500 

3.1 1 

1.93 

2.07 

7,2 

8.1 

Owen [21] 

Re =8000 

P 

12 

1.4 

1.43 

3.73 

4.0 


A, Pipe Expansion/Coaxial Jet: No Swirl 

The effects of various central jet/annular jet velocity ratios on the flow field can be reflected by 
the computed corner recirculation zone length X f which are listed in Table 1 The reattachment length 

for the coaxial jet with U a /U p = 1 is also shown. The wake created by the finite thickness of the pipe 

wall in coaxialjet flow extended to about one and a half chamber diameter downstream [Fig. 14(b)] , 
which produced extra longitudinal strain This can influence the length of the recirculation zone. As it 
can be seen from Table 1, an increase in the velocity ratio (U a /U p ) from 1 to 3 leads to a larger region 

of recirculation zone. Habib and Whitelaw [19] reported 10 percent increase in the length of the recir- 
culation zone for the higher velocity ratio (U„/U = 3) in their measurements. The k-e model prediction 

d. p 

14 
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Figure 13 Variation of centerline velocities. 


underestimates this increase (about 2 percent increase). For a similar experimental study by Johnson 
and Bennett [20] the calculated X f has a similar magnitude as the one calculated from Habib and 

Whitelaw's case, However, the measured X f of the latter is longer than the former. This may be attri- 
buted to the Reynolds number dependence as discussed in Section II, B, For the case of U a /Up =12, 

because of the large momentum difference, a large central recirculation zone was formed just downstream 
of the expansion. The strong shears created by the high velocity ratio of the annular and central jets 
largely enhance radial mixing thus reducing the region of the corner recirculation zone. This 
phenomenon is similar to that observed in the confined coaxial jet with strong swirl in the annular 
region (Section III.C). 

In Figure 13, the calculated centerline velocities for the confined coaxialjet with different ratios 
of the annular jet velocity to the central jet velocity are shown. The higher velocity ratio jets show more 
rapid mixing and decay compared to those for U a /Up = 1 in the axial direction. Also shown in Figure 

13 is the decay of centerline velocity of the pipe expansion flow In the region just downstream of the 
expansion, the decay of the coaxialjet with U a /Up = 1 is substantially faster than that of pipe expansion 

flow The experimental data of Reference 19 are also shown. The predictions are in fairly good agree- 
ment with the measurements. As suggested by Reference 19, the overprediction of mean velocity for 
unity velocity ratio and underprediction for a velocity ratio of three probably stem from the incorrect 
representation of the turbulent-diffusion process in the k-e model. Other possible discrepancies may be 
caused by the specification of the inlet values of k and e which have appreciable influence on the turbu- 
lent-diffusion process in the region downstream of the expansion. Figure 14 shows the normalized mean 
axial velocity profiles as a function of the radial location at different axial stations, X/D 0 , along the 
confined chamber for cases U a /Up = 1 and U a /Up = 3. At the plane just downstream of the expansion 
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(X/Do = 0. 1), the faster decay of the central jet for the velocity ratio U a /Up = 3 is clearly caused by 
the strong shear created by the velocity difference between Up and U a The mi nimum velocities in the 

wake region are due to the finite thickness of the inner pipe wall. The extent of the wake region is 
shorter for the case of higher velocity ratio (U a /Up = 3) because of the rapid mixing in the radial direc- 
tion. The predicted radial distributions of mean velocity compare fairly well with the measurements in 
the downstream region. However, the inaccurate velocity peak predicted at X/ D 0 = 1 43 indicates the 
deficiency of isotropic eddy viscosity hypothesis used in the k-e model. 

The longitudinal turbulence intensities normalized with Up are shown in Figure 15. The pre- 
dicted fluctuating quantities are derived from the isotropic turbulence assumption of k-e model, thus 

*»-(! ) m 


Ua/Up = 3 



Figure 16. Contours of kinetic energy. 
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Just downstream of the expansion, the maximum fluctuation occurs in the separated region. 
However, further downstream of the expansion, even after the reattachment point, the maximum fluc- 
tuation does not move close to the wall as would be expected. The pattern of the results can be further 
observed from the kinetic energy contours shown in Figure 16. These results show a distinct difference 
compared to the axi-symmetric pipe expansion flow (Figs. 6 and 7). The difference is likely attributed 
to the differences in the grometies, and that the confined coaxial jet has not reached the fully developed 
pipe flow at the exit plane of the calculation domain. The predicted axial fluctuations show correct 
qualitative trend for the confined coaxial jet without swirl, observed in experiments. The discrepancies 
in the upstream region predictions could well be due to the inadequate specification of inlet boundary 
conditions at the exit of central and annular jets. The sensitivities of the inlet boundary conditions on 
the development of the non-swirling coaxial jet will be addressed in the next section, 


B. Sensitivity to Inlet Boundary Conditions of Non-Swirling Flows 

In this section, the influences of inlet plane boundary conditions on calculated values of 
dependent variables are presented. The test case chosen is the confined coaxial jet with a sudden expan- 
sion with U a /Up = 3 of Habib and Whitelaw [191 In addition to the assumed shape of the radial 

profiles of the mean axial velocity, different levels of turbulence quantities such as k and e which are 
normally not known at the expansion plane will be tested to indicate the sensitivity of the inlet boundary 
conditions. 

The first case tested was the influence of the assumed shape of the velocity profile at the expan- 
sion plane. The entry profiles of k and e were calculated from equations (1) and (2). The calculations 
show that the specification of uniform velocities at the annular and central Jets lead to downstream values 
of the centerline velocity which are about 1 percent higher than those obtained with a 1/7-law velocity 
profile The calculations also show an insignificant increase in X r (~1 percent) for the fully-developed 

initial-velocity profiles as compared to the uniform-velocity profiles, The influence of the initial-velocity 
profiles can therefore be neglected. 

In the following calculations, the nrean-axial velocities are assumed uniform across the inlet plane. 
Usually, distributions of turbulence quantities such as k and e are not available from the experiments at 
the inlet plane, but are necessary for the calculation scheme as inlet conditions. Different initial turu- 
lent kinetic energy and energy dissipation rate would affect the eddy viscosity, i.e. the momentum trans- 
fer process. This may sometimes change the flow characteristics further downstream. Figure 17 presents 
results using three distributions of turbulent kinetic energy corresponding to (T) a uniform distribution 

with k/U^ n 2 = 0.003, (2) a distribution with 0.003 at the central nodes and 0.009 at the three near wall 

nodes, this profile was used to simulate a kinetic energy distribution in a fully developed pipe flow, 
and (3) a profile taken from Reference 31 with an empirical formulation for the fully developed pipe 
and annular flows 


k = 


(0.035 U^) 2 ^2 + 8 j , 0 < r < Rj 

(0.035 U in >2[2 + 4^y + 4^y] , R s <r<R 0 


(3) 
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Figure 17 Influence of initial profiles of kinetic energy on mean velocity 

The corresponding values of dissipation rate were determined from equation (2). The results 
show that the influence of the inlet profiles of kinetic energy is significant The influence upon the 
center-line mean velocity of an increase in the normalized kinetic energy from 0.003 to 0.009 at the 
three near-wall grid nodes is to shift the minimum velocity toward the flow entrance because of the large 
eddy viscosity assigned near the wall region, which enhances the spreading and decay of the central jet 
in the initial region The difference is also reflected in the regions away from the centerline and the 
entrance plane. The predicted recirculation zone length is appreciably shorter (10 percent) and the 
negative velocities attain lower values («7 percent) for the case of "clipped" kinetic energy profiles. The 
underpredictions of the maximum velocities after 1 D Q mentioned in Section III. A is much improved as 
a result of the new boundary conditions in k However, the inlet level of turbulence has significant 
influence on the kinetic energy distribution along the chamber Too much turbulence is introduced as 
a result of too large a kinetic energy in the initial region (Fig. 18) 

Finally, the importance of inlet e is considered here Estimation of e at the inlet plane usually 
is achieved by the mixing-length formulation [equation (2)] Another procedure of extracting e found 
in the literature was based on the isotropic eddy viscosity formulation 


e = -C 


dU 


M t i 3r 

p u v 01 
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X/Do 


Figure 18. Influence of initial profiles of k on centerline turbulent kinetic energy 

if u'v' and k were measured at the inlet plane, Equation (2) requires the specification of a characteristic 
dimension related to turbulent eddy size, Q' The inlet distribution of length scale has a large effect on 
the flow field in the upstream region where the diffusion term is important, The test of the inlet e level 
was performed by varying the e within the range 0.5 e r < e r < 2e r> where e r is the reference level of the 

E calculated from equation (2). As seen from Figure 19, the level of the e has appreciable influence on 
the calculations. As a matter of fact, the specification of appropriate initial turbulence level (k and e) 
plays a role of equal important to turbulence-model corrections as will be seen in Section III.C Further 
testing on the e level was carried out assigning the different distribution of 2' across the inlet plane 
The first case was to change the value of Q' in the annular exit using R Q instead of the annular width 

(R -R ) [22] The second case mvolved using an empirical correlation of the length scale distribution 
0 s 

across the pipe and annular region, the expression of 2’ is given by [23] 



The effects are shown in Figure 20 on the centerline-velocity decay and in Figure 21 on the e 
profile along the chamber The different specification of length scale distribution has a large effect in 
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Figure 21 Influence of initial profiles of e on the rate of dissipation 

the upstream region where the diffusion terms have a large effect on the mean velocity, it has no sig- 
nificant effect in the downstream region where the diffusion effects are minimal. If the swirl is intro- 
duced to the flow field, the turbulent diffusion process will be compounded by additional streamline 
curvature and extra strains, the calculations of the confined swirling flows are described in the next 
section 


C. Confined Swirling Flows and the Sensitivity to the Inlet Boundary Conditions 

The introduction of swirl into the flow creates much faster mixing, caused by radial pressure 
gradients and increase in the turbulence generation The flow field in a confined swirling jet is more 
difficult to predict due to the streamline curvature effect produced by the swirl and the strong aniso- 
tropy of the turbulence. In the study of confined coaxial swirlingjets, the swirl is introduced through 
the annular flow This study is essentially an extension of the confined coaxial jet and the geometrical 
configurations are the same as the ones used in Section III.A The azimuthal velocity profile of the 
annular flow can be a solid body rotation, a free vortex structure or just simply flat swirl velocity profile 
depending on the design of the swirl vanes used at the entrance of the expansion The overall swirl 
strength at the inlet plane of the expansion chamber is usually characterized by the swirl number S, 
which is a non-dimensional number representing axial flux of swirl momentum divided by axial flux of 
axial momentum, times the equivalent nozzle radius. In a coaxial swirlingjet, S can be written as 
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When swirl is introduced in the annular flow in the confined coaxial jet, significant radial and axial 
pressure gradients are set up near the central jet exit. At high degrees of swirl (S > 0.4 in confined 
coaxialjet), these pressure gradients would result in axial recirculation in the form of a central toroidal 
recirculation zone which is characterized by the presence of low tangential velocities, high turbulent 
intensities and large energy dissipation rates. The forming of the central toroidal recirculation zone 
depends on many factors, in addition to the swirl number, for example, size of enclosure chamber, 
geometry of the enclosure chamber, design of swirl vanes, and the inlet plane boundary conditions suck 
as particular velocity profiles and turbulence intensities. The more detailed descriptions of the current 
research activities in swirl flows can be found in Reference 24. 

Computations for the confined coaxialjet with swirl were initially made with two configurations, 
(i) Habib and Whitelaw [25] and (ii) Roback and Johnson [26] . These two configurations were basically 
extensions of the non-swirling jet cases of Habib and Whitelaw [19] and Johnson and Bennett [20] 

(Table 1). In Habib and Whitelaw's case, tangential entrance was used to generate a forced- vortex-type 
tangential velocity profile with S = 0.23. For the case of Roback and Johnson, a free-vortex type tan- 
gential velocity was generated by swirl vanes at the upstream of the expansion plane with S = 0.46. 
Because of the different swirl generator and the different swirl strength involved, the development of the 
two flow field shows drastic differences. Shown in Figure 22 is a comparison of the calculated mean 
axial velocities along the centerline with the experimental data. For the case of S = 0.46 a central 
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Figure 22. Variation of centerline velocityS = 0.23 and S = 0 . 46 . 
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recirculation zone appears in the initial region of the central jet. The present calculations with the grid 
23 x 24, which is close to the grid used by a similar study [14] , show very poor agreement with the 
measured data for the case S = 0.23. Further refinement in the grid system does not show any improve- 
ment. Some difficulties were encountered in a similar study by Sturgess [27] In his study, Sturgess 
indicated that fair comparisons could be obtained by careful placing of the gridlines, and the poor per- 
formance of the calculations resulted from the shortcomings of the physical turbulence model. Habib and 
Whitelaw [25] pointed out that the intermediate swirl number of 0.23 provides a flow which is probably 
more difficult to predict in that the near recirculating region is located away from solid surface and 
expansion plane. Different inlet conditions have been used as described in Section III.B trying to 
improve the calculations of centerline velocity distributions. With the higher level of kinetic energy pre- 
scribed for the inlet stream, which is shown to be very effective for the non-swirling case (Fig. 17), only 
in the very upstream region X/D Q < 0.5 the improvement is seen (about 5 percent); it has no significant 

effect in the downstream region (~ X/D Q >1) The serious discrepancies of the calculations for the 

intermediate swirl number case stem from the defects of the k-e model which is incapable of taking 
anisotropic and streamline curvature effects in the predictions. 

With a high swirl number, the recirculation region is closely tied to the expansion geometry and 
the downstream flow is more easily presented by the calculation methods. The subsequent calculations 
of the swirling confined coaxial jet with associated inlet boundary condition sensitivity test thus are pre- 
sented for the high swirl number case ( S s 0.46). Figure 23a shows the results using the distributions of 
kinetic energy corresponding to equation (1) and equation (3). The influence of the inlet stream turbu- 
lence on the centerline velocity distribution is not as significant as the non-swirling case (Fig. 17). The 
diffusion process induced by the strong swirl is more important than that involved by the inlet turbu- 
lence level. Thus the specification of the inlet tangential velocity profile is extremely important for the 
calculations of swirling flows. In fact, Ramos [5] showed that by using a Rankin-vortex type tangential 
velocity profile at the inlet plane of a swirling coaxial jet without expansion, k-e model gave drastically 
different predictions compared to the forced-vortex type tangential velocity specification at the inlet 
plane for both co- and counter-swirl flow conditions. A similar trend is found for the different specifica- 
tion of initial dissipation rate. In Figure 23b, the calculations were made with the inlet distributions of 
length scales estimated from equation (2) and equation (5) Also shown in Figure 23 are the predictions 
made with the modified k-e model based on the flux Richardson number correction of the e-equation 
[28] The effect of inlet e level is appreciable only in the initial region In a free-swirling j et study of 
Leschziner and Rodi [6] , the influence of the initial turbulence level was found to extend to the far 
downstream. In the confined-swirling flow, this influence is not significant in the downstream region 
In the upstream region, however, the specification of appropriate initial turbulence level (k and E) plays 
a role at least as important to turbulence-model corrections as seen in Figure 23, 


IV, COMPARISON OF THE PREDICTION OF PHOENICS AND TEACH CODES 


PHOENICS and TEACH codes use the same form of the time averaged equations governing the 
flow and the k-e turbulence model in which the differential equations for turbulence kinetic energy, k, 
and its dissipation rate, e, are employed. These equations are expressed in a general form and the appro- 
priate source term for each dependent variable is added [10] The equations are solved using finite 
difference techniques. 

Both the codes employ finite volume method of discretization in which integration is carried out 
over an elementary domain bounded by distinct lines of network. PHOENICS code used in the present 
computations employs an upwind difference scheme while TEACH employs a hybrid, central-upwind 
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scheme. Near wall treatment is through the wall function method. However, details of the design and 
execution of the solution algorithm may differ One wishes that the codes would predict correctly the 
physical problem at least qualitatively and would produce quantitative results as close to the measure- 
ments as possible. 

During the present investigation some discrepancies between predictions of the two codes have 
come to light. These are briefly described here to let the users of these codes know the deficiencies. 
Efforts should be made to improve the algorithms. The predictions compared here are for the flow in a 
sudden pipe expansion with an expansion ratio of 2, 

In a sudden pipe expansion, high shearing rates are generated in the separated layer (away from 
the wall). In this. region turbulence energy generation rate will be large due to the high shear, while 
dissipation rates will be low because they are away from the wall. The turbulence energy level thus 
increases rapidly to magnitudes many times larger than the regular pipe flow. The magnitude of the 
energy reduces further downstream and beyond reattachment point the flow develops to regular pipe 
flow. Thus, as discussed in Section ll,just downstream of the expansion, the kinetic energy maximum 
should occur in the separated layer, further downstream the maximum moves close to the wall. After 
the reattachments, the flow redevelops and attains the full developed profile. 

The kinetic energy contours predicted by the two codes are shown in Figure 24a, It is seen that 
the PHOENICS prediction produces kinetic energy contours in qualitative agreement with measurements, 
while the prediction of TEACH produces qualitatively different distribution of kinetic energy in some 
regions of the flow The reason for this is not known However, such contours have also been reported 
in the literature (see for example, Gosman, et al. [29] (Fig. 10a), TEACH prediction shows that the 
k-maximum moves .quickly toward the centerline of the pipe. Detailed examination of the algorithm is 
essential to explain this behavior of k and to incorporate modifications to rectify this problem. 

Another comparison is that of the important quantity, the size of the recirculation region. Both 
the codes predict the length of the reattachment more or less closely However, the locus of flow 
reversal (U = 0 contour) line is not predicted well. In addition, the PHOENICS code produces a long tail 
for the recirculation region. Physically such a long thin recirculation region cannot exist. If the tail is 
ignored the predicted length X r will be significantly smaller than the experimental value. Since 

PHOENICS code uses an upwind scheme, the false diffusion can influence the size of the recirculation 
region, It is not clear whether that is the reason for the observed results. The long tail region is also 
exhibited in the PHOENICS Demonstration Report by Qin [ 15] 

The implementation of the wall functions seems to be different in both the codes. As discussed 
in Section 2, the nodal point adjacent to the wall has the highest return flow velocity and at the next 
point the velocity drops down sharply This occurs in both codes. However, TEACH produces a 
smoothed velocity profile while PHOENICS produces a highly peaked profile as shown in Figure 24b. 

In TEACH code, special efforts are required to include blockages/irregular geometry etc. Hence, 
most of the computations of expansion flows, coaxialjets etc. using TEACH code, start the computation 
at the expansion plane. Thus, the important upstream effects on the flow as it arrives at the expansion 
are lost. Specifically, the kinetic energy profile at the expansion plane in both cases would be distinctly 
different. This has significant influence on the flow downstream. For example, the predicted reattach- 
ment length can differ by as much as one step height, i.e., is more than 10 percent, in the pipe expansion 
flow As would be expected, the computations start upstream and predict a larger X r The Xp thus 

predicted agrees more closely with the experimental measurements, since experimenters usually provide 
“enough” starting length. 
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(b) Flow reversal and velocity profile. 

Figure 24 Comparison of PHOENICS code and TEACH code, 





Thus, it is observed that the widely used PHOENICS and TEACH codes require some modifica- 
tions to produce predictions which would agree more closely with the experimental observations. 


V. CONCLUSIONS 


1 ) The severe underprediction of the length of the separated region in the case of the flow over 
a backward-facing step is related to the stronger streamline curvature in the initial region of the separated 
layer as compared to the axisymmetric flow 

2) In the case of the flow in a sudden pipe expansion, a larger inlet length scale produces a 
shorter reattachment length, as expected. 

3) The predicted reattachment length becomes independent of Re only for Re > 10^ 

4) At higher expansion ratios, the stabilizing effect of the wall in the symmetric expansion is 
significant. 

5 ) The finite thickness of the dividing lip at the expansion plane of the confined coaxial flows 
has significant influence on the flow development The large difference of the corner recirculation zone 
length between the coaxial jets with high velocity ratio (U a /Up = 3) and unity velocity ratio (U a /Up = 1) 
is underestimated by the k-e model prediction 

6) The calculations show that the inlet profiles of k and e have strong influence on the flow field 
of a confined jet without swirl. The effects of inlet boundary conditions extend to the further down- 
stream of the reattachment region. 

7 ) The influence of the inlet turbulence level in the confined coaxial swirling jet is less important 
compared to the non-swirling case. The influence is only appreciable in the initial region of the flow and 
has no significant effect in the downstream region, 

8) The k-e model failed to predict the momentum transport for swirling flows with intermediate 
swirl number Higher order turbulence models which take into account the anisotropic turbulent 
viscosity effects are necessary for the swirling flow field predictions, 

9) The k-e model prediction of PHOENICS and TEACH codes exhibit some qualitative discre- 
pancies in the predicted results when compared to the measurements. These codes may have to be modi- 
fied to improve the predictions. 


POSTSCRIPT 


Subsequent to the completion of this study, it was found that the qualitative difference on the 
predictions of the kinetic energy profiles discussed in Section IV was caused by the different averaging 
processes [c.f. Ref 10, p. 45] employed by the two computer codes. TEACH code uses arithmatic 
averaging while PHOENICS code utilizes harmonic averaging processes. Using an arithmatic averaging 
process, PHOENICS code produces kinetic energy contours and velocity profiles similar to those of the 
TEACH code. However, the long tail of the recirculation region observed when using PHOENICS code 
still persists. The tail region is probably caused by the different implementation of the wall functions. 
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According to Patankar [ 101, in the vicinity of wall region a rather fine grid, which was used in 
this study, would reduce the difference between the two averaging processes. The drastic difference of 
the kinetic energy contours produced by the two averaging processes near - the wall region (as discussed in 
Section IV) definitely needs further investigation. 
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